# Replication Material for "Social Capital, Institutional Rules 
#     and Constitutional Amendment Rates"
# William Blake, Joseph Cozza, Dave Armstrong and Amanda Friesen
#
# before running, execute setup.r first
#
# This file produces tables Table A8 and Figures 3 and A8. 

state_prop <- import("state_prop.dta")
state_rat <- import("state_rat.dta")

mod1n <- brm(bf(amendprop ~ vp_b + sc_b + lnwords + vp_w + sc_w +(1 | state),
                shape ~   sc_w +   (1  |state)) ,
             data=state_prop, 
             family=negbinomial(), 
             cores = 4, 
             backend = "cmdstanr", 
             refresh=0, 
             silent=2,
             seed=207, 
             save_pars = save_pars(all=TRUE))

mod2n <- brm(amendrat ~ vp2_b + sc_b + lnwords + sc_w + offset(log(amendprop)) + (1 |state),
             data=subset(state_rat, amendprop > 0), 
             family=negbinomial(), 
             #           control = list(adapt_delta=.9), 
             warmup = 2000, 
             iter=4000, 
             cores = 4, 
             backend = "cmdstanr",
             refresh=0,
             silent=2,
             seed=207, 
             save_pars = save_pars(all=TRUE))
mod1nb <- brm(bf(amendprop ~ vp_b + sc_b + lnwords + vp_w + sc_w +(1 | state),
                 shape ~   sc_w +   (1  |state)) ,
              data=subset(state_prop, ballotlimit2 >10), 
              family=negbinomial(), 
              cores = 4, 
              backend = "cmdstanr", 
              refresh=0, 
              silent=2, 
              seed=207, 
              save_pars = save_pars(all=TRUE))

mod2nb <- brm(amendrat ~ vp2_b + sc_b + lnwords + sc_w + offset(log(amendprop)) + (1 |state),
              data=subset(state_rat, amendprop > 0 & ballotlimit2 >10), 
              family=negbinomial(), 
              #           control = list(adapt_delta=.9), 
              warmup = 2000, 
              iter=4000, 
              cores = 4, 
              backend = "cmdstanr",
              refresh=0,
              silent=2, 
              seed=207, 
              save_pars = save_pars(all=TRUE))


sn1 <- summary(mod1n)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(obs = c(1,7,4:6, 2:3,8), 
         p = unname(c(sigf(mod1n)))) %>% 
  arrange(obs) %>% 
  select(-obs)

sn1b <- summary(mod1nb)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(obs = c(1,7,4:6, 2:3,8), 
         p = unname(c(sigf(mod1nb)))) %>% 
  arrange(obs) %>% 
  select(-obs)


sn2 <- summary(mod2n)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(obs = c(1,3:5, 2), 
         p = unname(c(sigf(mod2n)))) %>% 
  arrange(obs) %>% 
  select(-obs)

sn2b <- summary(mod2nb)$fixed %>% 
  select(c(1,3,4)) %>% 
  as_tibble(rownames="term") %>% 
  mutate(obs = c(1,3:5,2), 
         p = unname(c(sigf(mod2nb)))) %>% 
  arrange(obs) %>% 
  select(-obs)


s1r_c <- summary(mod1n)$random$state %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "State RE", term), 
         term = gsub("shape", "State RE (Shape)", term),
         p=unname(c(sigr(mod1n, "state"))))


s1r_cb <- summary(mod1nb)$random$state %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "State RE", term), 
         term = gsub("shape", "State RE (Shape)", term),
         p=unname(c(sigr(mod1n, "state"))))

s2r_c <- summary(mod2n)$random$state %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "State RE", term), 
         term = gsub("shape", "State RE (Shape)", term),
         p=unname(c(sigr(mod2n, "state"))))


s2r_cb <- summary(mod2nb)$random$state %>% select(c(1,3,4)) %>% as_tibble(rownames="term") %>% 
  mutate(term = gsub("Intercept", "State RE", term), 
         term = gsub("shape", "State RE (Shape)", term),
         p=unname(c(sigr(mod2nb, "state"))))

s1 <- bind_rows(sn1, s1r_c) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="vals") 

s2 <- bind_rows(sn2, s2r_c) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="vals") %>% 
  rename(vals2 = vals)

s1b <- bind_rows(sn1b, s1r_cb) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>%  
  pivot_longer(c(est, ci), names_to="param", values_to="vals")  %>% 
  rename(vals1b = vals)

s2b <- bind_rows(sn2b, s2r_cb) %>% 
  setNames(c("term", "estimate", "conf.low", "conf.high", "p")) %>% 
  mutate(sig = ifelse(p > .9, "*", ""), 
         est = sprintf("%.2f%s", estimate, sig), 
         ci = sprintf("(%.2f, %.2f)", conf.low, conf.high)) %>% 
  dplyr::select(term, est, ci) %>% 
  pivot_longer(c(est, ci), names_to="param", values_to="vals") %>% 
  rename(vals2b = vals)




outs <- full_join(s1, s2  %>% mutate(term = gsub("vp2_b", "vp_b",term))) %>% 
  full_join(s1b) %>% 
  full_join(s2b%>% mutate(term = gsub("vp2_b", "vp_b",term)))

outs$term[seq(2,20,by=2)] <- ""
outs <- outs %>% mutate(across(everything(), ~ifelse(is.na(.x), " ", .x)))

outs <-outs %>% 
  mutate(         term = gsub("shape_Intercept", "Intercept (overdispersion)", term), 
                  term = gsub("shape_sc_w", "Social Capital (within, overdispersion)", term),
                  term = gsub("vp_b", "Const. Rigidity (between)", term), 
                  term = gsub("vp_w", "Const. Rigidity (within)", term), 
                  term = gsub("sc_b", "Social Capital (between)", term), 
                  term = gsub("sc_w", "Social Capital (within)", term), 
                  term = gsub("tv", "Time Varying", term), 
                  term = gsub("lnwords", "log(# Words)", term)) 

outs <- outs %>% select(-param) %>% 
  bind_rows(data.frame(term = c("N", "N States"), 
                       vals = c("576", "48"), 
                       vals2 = c("452", "47"), 
                       vals1b = c("516", "43"), 
                       vals2b = c("405", "42"))) %>% 
  setNames(c("Term", "Proposal All Obs", "Ratification All Obs", 
             "Proposal No Limits", "Ratification No Limits"))


l1 <- loo(mod1n, moment_match=TRUE, reloo=TRUE)
l2 <- loo(mod2n, moment_match=TRUE)
l3 <- loo(mod1nb, moment_match=TRUE, reloo=TRUE)
l4 <- loo(mod2nb, moment_match=TRUE)

l_ests <- rbind(l1$estimates[3,], 
                l2$estimates[3,],
                l3$estimates[3,],
                l4$estimates[3,])


loos <- c("LOO IC (SE)", apply(l_ests, 1, function(x)sprintf("%.1f (%.1f)", x[1], x[2])))
names(loos) <- names(outs)
outs <- rbind(outs, loos)
r <- c("R(y, y-hat)", sprintf("%.2f", c(cor(fitted(mod1n)[,1], state_prop$amendprop),
                                        cor(fitted(mod2n)[,1], subset(state_rat, amendprop > 0)$amendrat),
                                        cor(fitted(mod1nb)[,1], subset(state_prop, ballotlimit2 >10)$amendprop),
                                        cor(fitted(mod2nb)[,1], subset(state_rat, amendprop > 0 & ballotlimit2 >10)$amendrat))))
outs <- rbind(outs, r)
outs

## Figure 3 (also Figure A5)

D1 <- D2 <- state_prop
s_sc <- sd(state_prop$sc_w, na.rm=TRUE)
D1$sc_w <- D1$sc_w - .5*s_sc
D2$sc_w <- D2$sc_w + .5*s_sc

fit1 <- fitted(mod1n, newdata=D1, summary=FALSE)
fit2 <- fitted(mod1n, newdata=D2, summary=FALSE)
diff1 <- fit2-fit1
diff1_ag <- t(diff1) %>% 
  as_tibble() %>% 
  mutate(state = state_prop$state) %>% 
  pivot_longer(-state, names_to = "iter", values_to = "vals") %>% 
  group_by(state, iter) %>% 
  summarise(m = sum(vals)) %>% 
  group_by(state) %>% 
  summarise(fit = mean(m), 
            lwr = unname(quantile(m, .05)), 
            upr = unname(quantile(m, .95)), 
            p = mean(m > 0)) %>% 
  arrange(fit)


state_prop$obs <- 1:nrow(state_prop)
D1 <- D2 <- left_join(state_rat, state_prop %>% select(state, year, obs))
#s_sc <- sd(tmp2$sc_w, na.rm=TRUE)

D1$sc_w <- D1$sc_w - .5*s_sc
D1$amendprop <- ceiling(colMeans(fit1))[D1$obs]
D2$sc_w <- D2$sc_w + .5*s_sc
D2$amendprop <- ceiling(colMeans(fit2))[D2$obs]

fit1 <- fitted(mod2n, newdata=D1, summary=FALSE, allow_new_levels=TRUE)
fit2 <- fitted(mod2n, newdata=D2, summary=FALSE, allow_new_levels=TRUE)
#diff2 <- fit2-fit1
f1 <- t(fit1) %>% 
  as_tibble() %>% 
  mutate(state = state_rat$state, 
         year = state_rat$year, 
         ap1 = D1$amendprop, 
         ap2 = D2$amendprop) %>% 
  pivot_longer(-c(state, year, ap1,ap2), names_to = "iter", values_to = "fit1")

f2 <- t(fit2) %>% 
  as_tibble() %>% 
  pivot_longer(everything(), names_to = "iter", values_to = "fit2")

diff2_ag <- f1 %>% 
  mutate(fit2 = f2$fit2) %>% 
  group_by(state, iter) %>% 
  summarise(across(c(fit1, fit2, ap1, ap2), sum)) %>% 
  group_by(state) %>% 
  summarise(fit = mean(fit2 - fit1), 
            lwr = unname(quantile(fit2 - fit1, .05)), 
            upr = unname(quantile(fit2 - fit1, .95)), 
            p = mean(fit2 - fit1 > 0), 
            ap1 = mean(ap1), 
            ap2 = mean(ap2)) %>% 
  mutate(apdiff = ap2-ap1) %>% 
  arrange(fit)

states_sf <- get_urbn_map("states", sf = TRUE)
tmp <- left_join(states_sf, diff2_ag %>% rename(state_name = state))
tmp <- tmp %>% filter(!is.na(fit))
col_blue <- brewer.pal(4, "Blues")
ggplot(tmp, aes(fill=log(fit))) + 
  geom_sf() + 
  theme_void() + 
  theme(legend.position="top", 
        plot.caption.position="plot", 
        plot.caption = element_text(hjust=0, size=12)) + 
  scale_fill_gradientn(colours=col_blue, 
                       breaks = log(c(1,3,9,27)), 
                       labels = c(1,3,9,27)) + 
  labs(fill="Predicted Additional\nAmendments Ratified", 
       caption = "See Table A8 (models 1-2) for model results")
ggsave("figures/fig3.png", height=8, widt=10, units="in", dpi=1000)
